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Abstract 

We propose a method for filling arbitrarily wide gaps in deterministic time series. Crucial to the 
method is the ability to apply Takens' theorem in order to reconstruct the dynamics underlying 
the time series. We introduce a functional to evaluate how compatible is a filling sequence of data 
with the reconstructed dynamics. An algorithm for minimizing the functional with a reasonable 
computational effort is then discussed. 

PACS numbers: 05.45.-a, 05.45.Ac, 05.45.Tp 
Keywords: time series analysis, chaos, filling gaps 



*Electronic address: francesco.paparella@unile.it 



1 



One problem faced by many practitioners in the applied sciences is the pres- 
ence of gaps (i.e. sequences of missing data) in observed time series, which 
makes hard or impossible any analysis. The problem is routinely solved by in- 
terpolation if the gap width is very short, but it becomes a formidable one if 
the gap width is larger than some time scale characterizing the predictability of 
the time series. 

If the physical system under study is described by a small set of coupled 
ordinary differential equations, then a theorem by Takens 3| suggests that 
from a single time series is possible to build-up a mathematical model whose dy- 
namics is diffeomorph to that of the system under examination. In this paper we 
leverage the dynamic reconstruction theorem of Takens for filling an arbitrarily 
wide gap in a time series. 

It is important to stress that the goal of the method is not that of recovering a 
good approximation to the lost data. Sensitive dependence on initial conditions, 
and imperfections of the reconstructed dynamics, make this goal a practical 
impossibility, except for some special cases, such as small gap width, or periodic 
dynamics. We rather aim at giving one or more surrogate data which can be 
considered compatible with the observed dynamics, in a sense which will be 
made rigorous in the following. 



I. INTRODUCTION 



We shall assume that an observable quantity s is a function of the state of a continuous- 
time, low-dimensional dynamical system, whose time evolution is confined on a strange 
attractor (that is, we explicitly discard transient behavior). Both the explicit form of the 
equations governing the dynamical system and the function which links its state to the signal 
s(t) may be unknown. We also assume that an instrument samples s(t) at regular intervals 
of length At, yelding an ordered set of N data 

Si = s((i-l)At), i = l,...,N. (1) 

If, for any cause, the instrument is unable to record the value of s for a number of times, 
there will be some invalid entries in the time series {sj}, for some values of the index i. 
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From the time series {sj} we reconstruct the underlying dynamics with the technique 
of delay coordinates. That is, we shall invoke Takens' theorem and claim that the 

m-dimensional vectors 



lie on a curve in M m which is diffeomorph to the curve followed in its (unknown) phase space 
by the state of the dynamical system which orginated the signal s(t). Here r is a positive 
integer, and i now runs only up to N = N — {m — l)r. Severals pitfalls have to be taken 
into account in order to choose the most appropriate values for m and r. Strong constraints 
also come from the length of the time series, compared to the characteristic time scales of 
the dynamical system, and from the amount of instrumental noise which affects the data. 
We shall not review these issues here, but address the reader to references 0,0,0]. 

We note that gaps (that is, invalid entries) in the time series {sj} do not prevent a 
successful reconstruction of a set 1Z = {xj} of state vectors, unless the total width of the 
gaps is comparable with N. We simply mark as "missing" any reconstructed vector Xj whose 
components are not all valid entries. 

If the valid vectors of 1Z sample well enough the underlying strange attractor embedded 
in M m , one may hope to find, by means of a suitable interpolation technique, a vector field 
F : U — > R m , such that within an open set U of M. m containing all the vectors Xj, the 
observed dynamics can be approximated by 



This very idea is at the base of several forecasting schemes, where one takes the last observed 
vector xjv as the initial condition for equation (J2J), and integrates it forward in time (see 
e.g., 0,3). 

The gap-filling problem was framed in terms of forecasts by Serre et. al. Their 
method, which amounts to a special form of the shooting algorithm for boundary value 
problems, is limited by the predictability properties of the dynamics, and cannot fill gaps of 
arbitrary width. 

The rest of this note is organized as follows: in section HU we cast the problem as a 
variational one, where a functional measures how well a candidate filling trajectory agrees 
with the vector field defining the observed dynamics. Then an algorithm is proposed for 
finding a filling trajectory. In section ITTT1 we give an example of what can be obtained with 
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this method. Finally, we discuss the algorithm and offer some speculations on future works 
in sec. IIV1 



II. A VARIATIONAL APPROACH 

The source of all difficulties of gap-filling comes from the following constraint: the inter- 
polating curve, which shall be as close as possible to a solution of (J2J, must start at the last 
valid vector before the gap and reach the first valid vector after the gap in a time T which 
is prescribed. 

To properly satisfy this constraint, we propose to frame the problem of filling gaps as 
a variational one. We are looking for a differentiable vector function £ : [0, T] — > U which 
minimizes the continuous functional 



Defining I = q — p, we have T = / At. If the curve coincided with the missing curve x(t) 
for t G [0, T], and F where a perfect reconstruction of the vector field governing the dynamics 
of the system, then the functional would reach its absolute minimum J = 0. The imperfect 
nature of F suggests that any curve which makes J small enough can be considered, on the 
basis of the available information, a surrogate of the true (missing) curve. In section Mil 
we shall offer a simple criterion to quantify how small is "small enough". For the moment 
we only care to remark that, even for a perfect reconstruction of the vector field, a curve £ 
making J abitrarily small, but not zero, need not to approximate x(£), in fact, the two curves 
may be quite different; however, such a curve £ is consistent with the dynamics prescribed 



Most methods for minimizing functionals focus on approaching as quickly as possible the 
local minimum closest to an initial guess. Because we may confidently assume that J has 
many relative minima far away from zero, we expect that these algorithms will fall on one 
of these uninteresting minima for most choices of the initial guess. Thus, our problem really 
reduces to that of finding a suitable initial guess. 




(3) 



with 



£(0) = x p ; f (T) 



by ©. 
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The complexity of the problem is greatly limited if we require that the set C = {y?}, 
whose Z + l elements sample the initial guess, has to be a subset of the set 7Z of reconstructed 
vectors. The index j = 0,1, ... ,1 does not necessarily follow the temporal order defined in 
1Z by the index i (cfr. eq. |H)), but we require that yo = x p and y; = x y . We shall denote 
with S(yj) the successor of the vector yj with respect to the temporal order in 7Z, and with 
P(yj) its predecessor. We propose that a good choice for the set C is one which makes as 
small as possible the following discretized form of the functional J: 



Of course, we shall restrict our choice of vectors to be included in C only to valid vectors of 
1Z having valid predecessor and successor. J could be zero if and only if £ cointained all 
the missing vectors in the correct order: C = {x p ,x p+ i, . . . ,x g }. The value of Jo increases 
every time that the order according to the j-index is different from the natural temporal 
order, that is, every time that S(yj) ^ Yj+i or P(yj) ^ Yj-i- If this happens we say that 
there is a jump in C between, respectively, the position j and j + 1, or j and j — 1. 

Although a set £ which performs many very small jumps may conceivably attain a very 
low value of Jo, there is an exceedingly small probability to find it within a finite dataset. An 
intuitive demonstration of this statement comes from the histogram in figure [TJ which shows 
the distribution of distancies between each reconstructed vector and its closest neighbor for 
the dataset discussed in section Mil as expected the frequency of closest neighbors quickly 
drops to zero for short distancies. Then our strategy will be that of looking for a set C which 
performs as few jumps as possible. 

Let us call orbit any sequence of valid vectors which does not jump. The first vector of 
an orbit shall have a valid predecessor, and the last a valid successor. Thus we define the 
predecessor of the orbit as the predecessor of its first vector and likewise the successor of the 
orbit as the successor of its last vector. We say that an orbit is consecutive to a point if its 
successor or its predecessor is the closest neighbor of the point. Two orbits are consecutive 
if the successor of one orbit is the closest neighbor of the first vector of the other orbit, or 
if the predecessor of one orbit is the closest neighbor of the last vector of the other orbit. 
Let us call branch a set made of consecutive orbits. Below we describe a simple algorithm 
to construct a set C by joining together one or more consecutive orbits. 

1. We follow forward in time the orbit consecutive to x p for I steps, or until it has a valid 
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successor. We store away the set of points made of x p followed by the points of this 
orbit as the I- jump forward branch. 

2. For each point of each (n — l)-jumps forward branch (where j = r,2r, . . . < If, r is 
an arbitrary stride, If + 1 is the number of points in the forward branch, and If < I,), 
we follow forward in time the orbit consecutive to for / — j steps, or until it has a 
valid successor. We store away all the points up to of the current forward branch 
followed by the points of the consecutive orbit as one of the n-jumps forward branches. 

3. We repeat step [2] for a fixed number rif of times. 

4. We follow backward in time the orbit consecutive to x g for / steps, or until it has a 
valid predecessor. We store away the points of the consecutive orbit followed by x g as 
the 1-jump backward branch. 

5. For each point zj of each (n — l)-jumps backward branch (where j — r, 2r, . . . < If,, 
r is an arbitrary stride, If, + 1 is the number of points in the backward branch, and 
h < I), we follow backward in time the orbit consecutive to Zj for j steps, or until 
it has a valid predecessor. We store away all the points of this orbit followed by all 
the points from zj to the end of the current backward branch as one of the n-jumps 
backward branches. 

6. We repeat step El a fixed number n& of times. 

7. For all possible pairs made by one forward branch and one backward branch we ex- 
amine synchronous pairs of points, that is a point y^. in the forward branch, and a 
point Zj b in the backward branch such that jf + h~ jb = I, where If, + 1 is the number 
of points in the backward branch. If they coincide, or one is the closest neighbour of 
the other, then we define C = {y , . . . , y jf ,z jb+1 , z h+1 }. 

The presence of jumps in C makes it unsuitable as a filling set for the gap. Furthermore, 
the discretized form (J3J of the functional Q is meaningful only if its arguments {y^} are a 
subset of the set of reconstructed vectors 71. Then we need a different discretization of Q 
which allows as argument any point of U. The simplest among many possibilities relies on 
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finite differences, leading to the following expression 



Ji ({Wi» 
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where the vectors Wj may or may not belong to 1Z, and j — 0,1, ... ,1. Here F(w J _ 1 / 2 ) is the 
vector field F evaluated at the midpoint between Wj_i and Wj. Ji is a function of m(l — 1) 
real variables (w = w p and wj + i = w q shall be kept fixed), which can be minimized with 
standard techniques, using £ as the initial guess. 



III. AN EXAMPLE 



In this section we show how the algorithm described above performs on a time series 

n 

generated by a chaotic attractor. We numerically integrate the Lorenz equations |H| with 
the usual parameters (a = 10, r = 28, b = 8/3). We sample the x— variable of the equations 
with an interval At = 0.02, collecting 5000 consecutive data points which are our time series. 
One thousand consecutive data points are then marked as "not-valid", thus inserting in the 
time series a gap with a width of l/5th of the series length, corresponding to a time T = 20. 
7 or this choice of parameters the Lorenz attractor has a positive Lyapunov exponent A ~ 0.9 
f|, setting the Lyapunov time scale at A" 1 ~ 1.1. We also find that the autocorrelation 
function of the time series drops to negligible values in about 3 time units. We conclude 
that T is well beyond any realistic predictability time for this time series. 

In the present example we selected the embedding delay r = 5 simply by visual inspection 
of the reconstructed trajectory, and we choose an embedding dimension m — 3. However, 
we checked that results are just as satisfactory up to (at least) embedding delay r = 15, and 
embedding dimension m = 6. 

We apply the algorithm with rij = 2 and n b = 0. The strides are r = 1 for the 2-jumps 
forward orbits and r = 100 for the 3-jumps forward orbits. This leads to 11001 forward 
orbits to be compared with 1 backward orbit, looking for synchronous pairs points which 
are neighbour of each other. We find two such pair of points, and the corresponding two 
initial guesses £1 and £2 are such that Jq(£i) = 1.21 and Jq(C2) = 1-04. 

The approximating vector field F is extremely simple, and its choice is dictated solely by 
ease of implementation. A comparison between different interpolating techniques is off the 
scope of this paper, and the interested reader may find further information in If Xj and 
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Xj are the vectors of 1Z, respectively, closest and second closest to Wj, then we define 

„, , x, — P(x;) + x, — P(Xo) , , 

F(wi-i/2) = ^ ■ (6) 

Using the definition © in we obtain Ji(£i) = 3.46 and ^(£2) — 3.17. In order to 
smooth-out the jumps in the filling sets, the function Ji is further decreased by iterating five 
times a steepest-descent line minimization (see, e.g., |U|) using £1 and £2 as initial guesses. 
This procedure yelds two sets of I + 1 points, .Mi and M.2 such that Ji(A'fi) = 1.54, and 
J\{M.-i) = 1.34. The corresponding time series are shown in figured The difference between 
the smoothed sets M. l and M.2 plotted in figure [2] and the sets with jumps £1 and £ 2 would 
be barely noticeable on the scale of the plot. 

The effect of the smoothing may be appreciated by looking at figure El which shows the 
region across the jump between two consecutive orbits of L\. The non-smoothed filling set 
(dashed line) abruptly jumps from one orbit to the other, but the smoothed trajectory (thick 
solid line) singled out by the points of M.\ gently moves between them. 

No attempt has been made to approach as closely as possible the local minimum of J\. 
In fact, we verified that for orbits in 71 having the same length as the interpolating sets A4% 
and A4 2l J\ ranges (roughly) between 1 and 9. This is a measure of the accuracy with wich 
the field F approximates the true dynamics of the observed system, and there is no point in 
looking for an interpolating set having a value of J\ below this range. 



IV. DISCUSSION AND CONCLUSIONS. 



In this note we have described an algorithm which fills an arbitrarily wide gap in a time 
series, provided that the dynamic reconstruction method of Takens is applicable. The goal 
is to provide a filling signal which is consistent with the observed dynamics, in the sense 
that, in the reconstructed phase space, the vector tangent to the filling curve should be close 
to the vector field modeling the observed dynamics. This request is cast as a variational 
problem, defined by the functional (j3J. The acceptable degree of closeness is determined by 
the level of accuracy of the reconstruction, which we quantify by computing the discretized 
form (jnj) of the functional for orbits having the same length of the gap. 

Obviously, if the time series has more than one gap, our method can be applied to all the 
gaps, indipendently from each other. 
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A second novel idea, that greatly simplyfies the problem, is that of building a rough 
solution by stitching together pieces of the observed dataset. The actual solution, which will 
not be an exact copy of anything present in the observed dataset, is obtained by refining 
this first draft. We have illustrated a basic algorithm that embodies this idea, although 
no attempt has been made at making it computationally optimal. In particular, with the 
algorithm in its present form, many of the forward and backward branches will be partial 
copies of each other, because nothing forbids two distinct branches to jump on the same 
orbit. This leaves some room for improvement, because the effectiveness of the method 
relies on a substantial amount of the set 7Z of reconstructed vectors to be explored by a 
relatively limited number of branches. In a forthcoming, enhanced version of the algorithm 
some kind of tagging mechanism shall be incorporated in order to produce non-overlapping 
hierarchies of forward and backward branches. 

We observe that this algorithm does not give a guarantee of success: it is perfectly 
possible that no point of the forward branches is the closest neighbor of (or coincides with) 
a synchronous point of the backward branches. In this case the obvious attempt is to 
deepen the hierarchy of the branches, as much as it is computationally feasible. Or, one 
may relax the request that branches may jump only between closest neighbors, and accept 
jumps between second or third neighbors as well. As a last resort, one may stitch any pair of 
forward and backward branches at their closest synchronous points, hoping that the resulting 
jump could later be smoothed satisfactorily by minimizing the function (JHJ). However, when 
facing a failure of the algorithm, we believe that first should be questioned the goodness 
and appropriateness of the dynamic reconstruction. The presence of too many gaps, the 
shortness of the time series, or measurement inaccuracies may make the gap-filling problem 
an insoluble one. We speculate that the ability of filling gaps with relative ease is a a way 
to test the goodness of a dynamic reconstruction. 

The ease with which a gap may be filled, as a function of his width, is a problem deserving 
further work. For the moment we simply recall that if a set of initial conditions of non-zero 
measure is evolved in time according to (j2J, eventually we expect it to spread everywhere 
on the attractor (here the measure is the physical measure /i of the attractor cfr. ref. (l^). 
More rigorously, if <f> t is the flow associated to (pj, and if it is a mixing transformation, 
then, for any pair of sets A, B of non-zero measure, lim^^ u{<f>tA PI B) — /z(A)/i(£>). The 
dispersion of a set of initial conditions is further discussed in where, for example, it is 
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shown that the essential diameter of a set of initial conditions cannot decrease in time, after 
an initial transient of finite length. 

This leads to the idea that wide gaps should be easier to fill than not-so-wide ones, 
because forward and backward branches have explored larger portions of the attractor, and 
so there is a greater chance to find synchronous points where they can be joined together. 
As a first step toward the verification of this hypotesis, we computed the average minimum 
distance between synchronous points of the 1-jump forward and backward branch as the 
gap moves along the dataset, for several gap widths. We used the dataset discussed in sec. 
IIIII and a ten times longer extension of it. The results, plotted in figure IH show that the 
average separation of the branches initially increases as the gap widens, but then it reaches 
a well-defined maximum and, from there on, decreases as the gap width is further increased. 

We close by mentioning that the dynamic reconstruction technique has been successfully 
applied even to stationary stochastic time series, to generate surrogate data with the same 
statistics of the observed ones This fact, and the hypothesis that ergodicity (or the 

stronger requirement of being mixing), rather than determinism, is the crucial property that 
allows for filling gaps, suggest that some modified version of our method should be able to 
fill gaps in a large class of stochastic time series. 
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Figure Captions 

Figure [T] Distibution of the distances between each reconstructed vector and its closest 
neighbor for the dataset discussed in sec. Mil 

Figure [2] Panel A) shows portion of the time series discussed in sec. Mil The blackened 
line was removed and the resulting gap was filled by applying the algorithm described in 
section HU The blackened lines in panels B) and C) are two different fillings. 

Figure [3] Plot of the first two components of the reconstructed vectors of: A4y (thick solid 
line with crosses); C\ (thick dashed line with asterisks); one orbit of C\ and its successors 
(thin line with open circles); the orbit consecutive to it and its predecessors (thin line with 
open squares). To illustrate the smoothing effect of minimizing functional (jSJ, we only plot 
a very small portion of these sets in the vicinity of the jump between the consecutive orbits. 

Figure |4] Average minimum distance of sinchronous points in the 1-jump forward and 
backward branch as a function of gap width. The dashed line refers to the dataset discussed 
in sec. MIL the solid line refers to a ten times longer extension of that dataset. The vertical 
line marks the Lyapunov time A^ 1 ~ 1.1. 
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Distance of the Closest Neighbor 



Figure 1: Distribution of the distances between each reconstructed vector and its closest neighbor 
for the dataset discussed in sec. II I II 
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Figure 2: Panel A) shows portion of the time series discussed in sec. IT77I The blackened line was 
removed and the resulting gap was filled by applying the algorithm described in section [HI The 
blackened lines in panels B) and C) are two different fillings. 
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First Delay Coordinate 



Figure 3: Plot of the first two components of the reconstructed vectors of: M\ (thick solid line 
with crosses); C\ (thick dashed line with asterisks); one orbit of L\ and its successors (thin line 
with open circles); the orbit consecutive to it and its predecessors (thin line with open squares). 
To illustrate the smoothing effect of minimizing functional (JHJ, we only plot a very small portion 
of these sets in the vicinity of the jump between the consecutive orbits. 
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Figure 4: Average minimum distance of sinchronous points in the 1-jump forward and backward 
branch as a function of gap width. The dashed line refers to the dataset discussed in sec. IIH[ 
the solid line refers to a ten times longer extension of that dataset. The vertical line marks the 
Lyapunov time A -1 ~ 1.1. 
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